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ABSTRACT 


Baroclinic instability in the ocean is a primary cause of mesoscale eddies, which 
are pockets of water in the scale of 100km that have different density, thennal, and 
rotational characteristics than their surroundings. First observed in the early 1900s, eddies 
are thought to be a predominant reason for the heat flux between the equator and the 
poles in both the ocean and the atmosphere. In attempt to understand this process better, 
this study uses a series of numerical simulations perfonned on high performance 
computing systems. The calculations are based on the Massachusetts Institute of 
Technology general circulation model, which is used to compare lateral heat transport 
between different simulations. The specific objectives of this project include 

i) Comparison the direct and remote interactions of shear with topographic slope. 
The direct scenario is one in which the shear extends throughout the entire ocean depth 
and is therefore in direct contact with the sea floor, whereas in the remote scenario there 
is a spatial separation between the shear in the upper half of the basin and the bottom 
topography, 

ii) Analysis of the system response to changes in the zonal and meridional 
seafloor slope, and 

iii) Investigation of the effect of orientation changes in the mean large-scale 
current on cross-flow fluxes. 

The lateral heat transport and diffusivity of these simulations are then compared 
to our analytic model, known as Growth Rate Balance, which is based on the balance 
between growth rate (primary) instabilities deduced from linear theory and numerically 
generated secondary instabilities. 
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I. 


INTRODUCTION 


Modern science has both the distinction of advancing at a very fast rate due to 
discoveries in technology and mathematics, and the ignominy of having advanced so 
quickly that society has become unaware of the many natural phenomena that still remain 
unexplained. This thesis explores one such topic whose effect on mixing has eluded 
scientists from multiple disciplines and continues to be an area of debate and research. 
The baroclinic instability (BI), which evolves into turbulent flows that transport heat and 
momentum, is a predominant source of mesoscale variability in both the ocean and in the 
atmosphere. The nonlinear interaction of flow components leading to growing 
disturbances creates large-scale flows that can be easily identified in nature through in- 
situ and remote observations (Kamenkovich et al. 1986). These disturbances break off 
into energetic, relatively long-lasting circular features known as eddies. 

A. BACKGROUND 

Eddies are everywhere! In the atmosphere, they can be seen in various forms, 
such as that of mid-latitude extra tropical cyclones. Mesoscale eddies are features that 
dominate the ocean basins, especially around strong oceanic current systems like the Gulf 
Stream or the Antarctic Circumpolar Current (ACC), and are scientifically important for 
many reasons: For one, eddies are believed to be the source of more energy than any 
other ocean motion. Accordingly, they also transport tracers, chemicals, and nutrients, 
making them a cross-disciplinary topic of study (Robinson 1983). One of the greatest 
aspects of eddies is the transfer of heat from the tropics to the poles; thus eddy heat 
transport is believed to be a major contributor to the balance of the global heat budget in 
both the atmosphere and the ocean. Without large-scale eddy transport in the atmosphere, 
for example, the mean meridional circulation would be too balanced to explain the actual 
bulk energy transport (Green 1970). Figure 1 illustrates a satellite-derived image of 
eddies in the Gulf of Mexico and in the Gulf Stream. 

Ocean eddies have spatial scales ranging anywhere from thousands to hundreds of 
thousands of meters, and time periods varying from weeks to years (Robinson 1983). 
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Because of this, parameterizations of eddies in climate models remain a major source of 
uncertainty which can lead to considerable inaccuracies in model predictions (Gent and 
McWilliams 1990). This feature is rather unfortunate, considering the importance of 
eddies to global heat budgets and in enhancing model accuracy in the long-range 
forecasts. For example, our inability to properly represent eddies serves as a potential 
cause of great uncertainty in ice melting in Arctic Ocean models (Maslowski and 
Lipscomb 2003). Unphysical eddy closure models can parameterize heat flux to the south 
for places like the Canadian Basin where northward (positive) heat flux is observed in 
nature (Steiner et al. 2004). The inadequate parameterization of eddies diminishes 
accuracy and customer confidence. 



NASA/Goddard Space Flight Center Scientific Visualization Studio http://svs.gsfc.nasa.gov/3827 


Figure 1. Perpetual Ocean Image. Source: NASA Goddard Space 

Flight Center (2011). 

Traditionally, studies of ocean eddies have been more theoretical than 
observational. Pre-satellite eddy observations were few and far between. Eddies were 
first documented by Ben Franklin’s grandnephew, Jonathan Williams, who noted a 
“warm core ring” while monitoring sea surface temperatures and velocities in the Gulf 
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Stream (Robinson 1983, Williams 1793). Observations for these “rings” was noted again 
in the 1930s by P.E. Church who discovered warm core eddies while examining ship 
thennograph records from Gulf Stream transits (Church 1932, 1937). Iselin (1936) 
followed this up with deep-water temperature and salinity samples taken from The 
Atlantis, from which he concluded these rings to be a pennanent feature. By the 1970s, 
teams of scientists from around the world teamed up on surveys such as MODE 1 and 
POLYMODE, which was the largest joint U.S.-U.S.S.R. experiment of its time 
(Robinson 1983). Now, with the use of satellites, eddies can be observed and tracked 
around the globe on a daily basis. Additionally, high perfonnance computer systems and 
GCMs (General Circulation Models) can realistically model the ocean in three 
dimensions with increasing speed and resolution. Although we may not know as much 
about mesoscale variability as we would like, technology has given us the advanced tools 
necessary that we may one day be able to solve many of Earth’s mysteries (Knauss 
2000 ). 

B. TOPOGRAPHIC EFFECTS AND SHEAR 

Of the many topics that require further investigation with regard to ocean eddy 
dynamics, the effects of topography on synoptic eddies is rapidly gaining attention of the 
oceanographic community. There are many ways that bottom topography can impact the 
physical ocean processes such as sea floor roughness, smooth versus rocky geologic 
features, and sea floor gradient. Observations indicate that roughness on sloped 
topography, even in the deep ocean abyss, will enhance vertical mixing and turbulent 
diffusivity (Polzin et al. 1997, Dewar 1998). This creates an issue with models that have a 
bathymetric resolution too course to represent roughness, which may cause inaccurate 
estimates of diffusivity and other flow properties, especially near the sea floor. This 
inaccuracy can be especially prevalent when attempting to model processes like internal 
waves or eddies (Goff and Arbic 2010). Thus, there is a need for good synthetic 
parameterization of this roughness, particularly with lower-resolution models. 


1 Mid-Ocean Dynamics Experiment. 
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Bottom relief, or slope, is known as a generating mechanism for synoptic 
variability. There are various ways slope can affect large- and meso-scale patterns. For 
example, topographic Rossby waves are a low-frequency phenomenon that can disappear 
when the wave is no longer in contact with the sea floor. Additionally, flow over bottom 
irregularities, such as sea mounts and guyots, can create instabilities which, due to the 
conservation of potential vorticity, can create cold core eddies above them with warm 
core eddies downstream (Kamenkovich et al., 1986). Topographic stress generated on 
these eddies can create secondary circulations which may increase upwelling and sea 
surface height (Holloway 1987). Sutyrin and Grimshaw (2005, 2010) considered the 
effects of sloping topography by applying surface-intensified circular vortices on a P- 
plane in a two-layer model with reduced gravity approximation. In their experiment they 
specifically looked at frictional effects on deep-ocean flows evaluating the topographic 
orientation effects using the p-drift with either along-slope or cross-slope components 
(Sutyrin and Grimshaw 2005). Our knowledge of the way various slopes affects eddy 
trajectory, for example, and how it changes the structure of the eddy as it forms, is still 
developing (Robinson 1983). 

A third type of effect caused by bottom relief is the impact it has on the stability 
of the current itself. As we will discuss in the next chapter, zonal currents can be 
maintained in the ocean without external forces. In a stratified system with a flat (uniform 
depth) seafloor and zonal current, shear that remains vertically unifonn can be 
completely stable with the only varying parameter being the Coriolis force: a fictitious 
force that changes with latitude as measured by the aptly-named P-effect, (1), 
(Kamenkovich et al. 1986). 



( 1 ) 


Keeping all the same background conditions, the addition of a bottom slope can 
create instability (Charney and Flierl 1981). Green (1970) also noted the importance of 
shear in meridional heat flux and that constant shear in the troposphere would create an 
entropy flux that was independent of height. In this, he suggested that when there is 

inconsistent shear in the vertical, the heat flux can change sign. If this holds true in the 
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ocean, we could expect that negative values of heat transfer exist below layers of current- 
induced shear. The remote case in our experiment is an example of inconsistent shear. 
This invites the question of how the slope of the topography in a motionless lower layer 
will impact the heat flux in the upper remote layer, and if there will be a conservation 
effect that will shift the sign of heat flux at particular depths, particularly below the level 
of no shear. 

In this study, we attempt to isolate the effect of slope on eddy dynamics by using 
high-resolution numerical simulations in a hydrostatically balanced, eddy-resolving 
model. We analyze the equilibrium dynamics of BI in an idealized environment and 
compare meridional diffusivity of a flat bottom case with the diffusivity over various 
slopes. Meridional slope should have a de-stabilizing effect on shear (Chen and 
Kamenkovich 2013) which should affect the eddy-transfer process. It will be shown that 
where slope exists, there is an increase in bottom interaction, and a decrease in stability, 
both of which enhance eddy formation. The question we look to answer is this: How will 
increased positive and negative (north and south) slopes affect large-scale eddy-induced 
transfer? 

C. GROWTH RATE BALANCE MODEL 

One of the great difficulties with our understanding of eddy dynamics is the non¬ 
linearity of BI. In recent decades, many attempts have been made to explain mesoscale 
variability in tenns of analytical models. Because mesoscale variability is non-linear, full- 
fledged models would be computationally prohibitive and mathematically difficult 
(McWilliams 2011). Several researchers, such as Thompson (2010) and Visbeck et al. 
(1997), have come up with different unique ways to simplify this reasoning using analytical 
models. One of these theoretical models, which was developed in 2014, is known as the 
growth rate balance model (GRB). This model is based on the assumed balance between 
primary (/./) and secondary (fy) linear modes of instability, which are linked with an 
empirical constant (C). This paper will compare our 3-D model results with those of the GRB 
model (Radko et al. 2014). 
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D. ORGANIZATION 


The following chapters of this thesis will contain the model description as it is 
used in MITgcm. This includes the stability theory that explains how the model creates 
eddies. We will then explain the initial conditions and boundary conditions, rationale for 
the conditions that were chosen, and differences between the types of model runs that 
were conducted in this experiment. Afterward, the diagnostic tools will be discussed 
before moving on to Chapter III, numeric results, which will include data, charts and 
snapshots that best represent the data and new discoveries. Following that will be our 
comparison with GRB theory results in Chapter IV, discussion and conclusions in 
Chapter V then recommendations for future research in this topic in Chapter VI. 
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II. MODEL DESCRIPTION 


A. FORMULATION 

An important step in modeling ocean eddies is to understand how BI is created 
within the model. This chapter gives us a description of the model dynamics, explains the 
assumptions, and outlines rationale behind input parameters. 


1. Hydrostatic Equations 

In this project, we will look at BI in terms of a hydrostatic, continuously stratified 
model with a uniform horizontal velocity decreasing with depth. Of the three options, 
hydrostatic, quasi-hydrostatic and non-hydrostatic, the former was chosen for the numeric 
representation. Although the non-hydrostatic model is the preferred choice because it is 
the most general, it is also the most complex and processor intensive for a numeric 
model. The hydrostatic model can be as effective as either of its counterparts in large- 
scale scenarios (Marshall et al. 1997). By scaling the dimensions of the model, a 
determination can be made that hydrostatic approximation is valid. Our length scales are 
2000 km, thus L = 2xl0 6 m. The height of our basin is 3.3 km, making H ~ 3xl0 3 m, our 
average meridional velocity is 0.075 m-s' 1 in the remote case and 0.15 m-s' 1 in the direct 
which we can approximate U ~ 0.1 m-s' 1 . Using a standard sea-water potential density of 
I x 10 kg-m' , which means that—with a temperature difference from surface to bottom 
of five degrees—we can assume a change in density over depth (based on temperature) of 
1 kg-m' . Therefore, we can calculate our buoyancy frequency using Equation (2) below 
to be N ~ 1.7xl0' 3 s' 1 . 


N 2 = 


S_ dp 

p dz 


( 2 ) 


9 . 

Here, g is the standard gravitational constant of 9.8 m-s' and p is water density, which 
will be described in Section B. In the direct case, the velocity structure of the basic state 
is as follows with v representing meridional velocity and u representing zonal velocity 
with the subscripts indicating their vertical location in the water column. 
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( 3 ) 


V — 0 9 14 surface — U x9 U bottom — 0, 
d I t (tl surface- U bottom ) 

dz ~Z 


In the remote case, however, these expressions change to 


— — z 

Usurface = W m , U f . (/? = —) = 0, 


Oil U surface Ur 

dz 1 


(4) 


In both cases, remote and direct, the velocity fields are related to the temperature patterns 
through the thermal wind balance: 


dT _ f du 
dy ga dz 


(5) 


Initial boundary conditions were adopted based on Radko et al. (2014) with 
changes made to reflect a deeper open-ocean seafloor that would better represent 
conditions in the mid-ocean basins. This is particularly pertinent to the Southern Ocean, 
which reaches depths over 5 km. Because the model setup had already been tested, this 
became the basis for this experiment. This experiment includes /? of lxlO' 11 and standard 
Coriolis parameter of 1x10 4 nf’-s . Density stratification in our model ocean basin was 
created by assuming uniform salinity of 35 psu, thereby making density solely a function 
of temperature. Velocity is induced via Equations (3) and (4). In order to create both a 
horizontal and vertical stratification, the equation for thermal wind (5) is used in 
conjunction with shear gradient to generate a temperature gradient in the model 
(Kamenkovich et al. 2009). This temperature gradient corresponds to an approximate 
10°C overall variation in the vertical and 14°C in the horizontal. A randomly generated 
temperature variation of magnitude 0.10°C was induced throughout the model, which 
disrupts the initial stratification enough to facilitate instability once it is acted on by 
shear. 
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2 . 


Shear 


One of the main drivers of BI is the shear. Vertical shear provides a key 
component into destabilizing the stratification and thereby inducing instability. For this 
study we chose a magnitude of shear that is observationally reasonable capable of 
generating BI. By testing standard surface current velocities, u surface , with a gradient of 

decreasing velocity with depth,--, a standard gradient was determined that is used in 

dz 

all our shear and depth scenarios. In this research, two shear configurations and two 
depths are compared. The first configuration is the direct case where the shear was 
induced throughout the water column, reducing linearly in intensity with depth until 
u bottom = 0. The second configuration is the remote case, which is where the shear was 

applied to the upper half of the basin, but the lower half was given motionless initial 
conditions. No meridional flow is initially induced and therefore any flow in the y- 
direction is caused by dynamic effects. In doing this, a comparison can be made between 
both the direct and remote shear scenarios and between the different meridional slopes. 
By measuring the eddy heat flux, we can determine if there is a remote effect on the 
topography on the background current, and the role that bottom slope plays on ocean heat 
transport. 

The remote case scenario was conducted at two different depths: 1100 m and 
3300 m. In order to compare the shear effects between scenarios at different depths, we 
required that the shear gradient remain the same instead of the mean or maximum 
velocities of zonal flow. Therefore, for the 3300m case, the surface velocity measured 
0.15 m-s' 1 and reached zero velocity at a depth (z) of 1550 m whereas the 1100 m case 
measured 0.05 m-s' 1 reaches zero velocity at 550 m. This gives give us a shear gradient 
that is equivalent for both depths, despite a greater area covered by shear in the 3300-m 
case. Thus, our comparisons are based strictly on a comparable shear gradient rather than 
separate shear gradients with the same surface velocity as illustrated by 

p-jji 

—(h 0 = 1100m) = 0.909le ^ 1 =—(h 0 = 3300/n). (6) 

dz dz 
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3. 


Zonal versus Non-Zonal Currents: Orientation 


Because of the Coriolis force, a current in the zonal direction is the only flow that 
can be maintained without the effects of outside forces (Kamenkovich et al. 1986). In this 
experiment, the flow direction is turned at 30° intervals over a flat bottom basin, and the 
impact of the rotation on the cross-current thermal flux (CFF) is analyzed. The magnitude 

of horizontal velocity is U = yju 2 + V 2 , where U represents along-flow current and V 
represents the across-flow current. The term “flow” in this context describes the initial 
background current in the direction of orientation. Lower case x and y represent the 
longitudinal and latitudinal axes, which is standard in our zonal flow simulations, but in 
the orientation cases we use X and Y to indicate the along-flow and cross-flow axes 
respectively. Because the dimensions of the basin remain the same, and the basin-relative 
flow remains the same, the value of any diagnostic in Y, although no longer meridional, is 
our desired quantity. The thennal gradient also changes with each rotation, and therefore, 
it is the (3-effect that will provide the physical difference between orientation runs. 

It should be emphasized that, in the orientation cases the time mean flow is not 
strictly oriented in the X-direction. Therefore, the mean temperature flux in Y includes a 
cross-flow component associated with the advection of heat by time-mean velocity, 
which must be accounted for in order to determine CFF. This is done in the post¬ 
processing analysis as follows: 

rV' = TV-TV. (7) 

B. INITIAL AND BOUNDARY CONDITIONS 

1. Size and Parameters 

In detennining the size of the basin, a few principles of eddy dynamics were 
implemented. Since mesoscale eddies are in the scale of tens to hundreds of kilometers 
wide, the horizontal dimensions of the box had to be large enough to resolve 100 km 
features while the resolution also had to be fine enough to encompass the smaller eddies. 
In particular, baroclinic radius of deformation ( R d ) must be resolved in order to ensure 
that baroclinic instability is properly represented. Inherently, the Brunt-Vaisala buoyancy 
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frequency (AO determines radius of deformation. To simplify the equation for detennining 
density stratification, salinity was kept constant at 35 psu, which means that density(p) 
becomes a function of temperature: p = p(T). Because of this, there is a direct relationship 
between our temperature gradient and pressure gradient, which can be used to estimate 
buoyancy. The expressions below use the equation of state to replace density with 
temperature and solve for R d . Note that a is the thermal expansion coefficient. 


R d = 


Nh_ 

f 


dS = Q,^-^ = a 
Po 


[T~T„) 




dp 

dz 


dT_ 

dz 


( 8 ) 

( 9 ) 


l~g d P 
V Po dz 



( 10 ) 


This equation gives a calculated R d of approximately 2.4x 10 4 m, or 24 km for the 

1100-meter runs and 4.20x 10 4 m (42 km) for the 3300-meter case. This is consistent with 
the size of mesoscale eddies. This estimate led to the optimal choice for basin size for this 
experiment of 2000 km by 2000 km, as the basin had to be wide enough to a contain a 
large number of eddies. Additionally, the horizontal resolution is also a factor since too 
course a resolution will not be able to simulate eddies, and too fine a resolution will cause 
model run times to increase dramatically. Two-kilometer resolution was chosen as it is 
1/10 th the scale of concern of R d for the 1100-meter case. 

2. Direct and Remote Case 

As mentioned in the introduction, two major flow regimes were introduced into 
the original research plan, which have been aptly named the direct and remote cases. The 
direct case is perhaps the more straightforward of the two as the shear is exerted 
throughout the water column by introducing a surface current that decreases in magnitude 
linearly with depth all the way to the seafloor. Thus, the shear is present throughout the 

water column and is only inhibited at the sea floor. In the remote case, the shear is 
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present only in the upper half of the water column, and the velocity becomes zero at half 
depth. Thus, the lower half of the water column is completely motionless initially. Figure 
2 is a visualization of the remote and direct shear scenarios. A 3-D visualization can be 
seen in Appendix A (Figure 18). 


DIRECT 


REMOTE 



Figure 2. Diagram of the Direct (left) and Remote (right) Shear Scenarios 


This motionless region is not a “lower layer” and should not be confused with a 
two-layer model because the stratification in the upper and lower half is continuous. The 
lower half is merely a region of no initial zonal current, and therefore, any flow that 
occurs in the lower layer is due to baroclinic instability effects induced in the upper layer. 
Also, the magnitude of the shear is the same across the direct and remote cases is the 
same. That is the linear decrease in velocity is the same. This also means that the velocity 
on the surface is twice the value in the direct case. This allows us to maintain the same 
shear magnitude despite having larger surface speeds with the direct flow: 


r)ii 

* ir (h f) = 1100m) « 0.055m.y ', ^ d,r (h n = 1100m) « 0.026mA 


- 1/2 


= 3300m) * 0.15 ms~\j^(h 0 = 3300m) * 0.078mA 1 . 


( 11 ) 


- 1/2 
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3. 


Zonal Case 


The previous model compared several cases of various meridional, or north— 
south slopes. In an effort to mimic the somewhat of a mid-Atlantic ridge, a series of runs 
were performed to determine the effects of topographic slope in a zonal, or east—west, 
direction, as depicted by Figure 3. In these scenarios, the depth of the water remained 
constant in y, but changed over x with an apex at length x/2; thus, the flow is 
perpendicular to the bottom gradient. This also means that potential vorticity of the mean 
flow in the direct case will obviously be affected by the height of the water column as 
seen by Equation (12). However, the impact of the bottom topography in the remote case 
will be less certain. Additionally, because the zonal slope can deflect the initially zonal 
flow, the background meridional flow will have to be taken into account when evaluating 
the eddy-induced transport. Equations for the potential vorticity and zonal slope are 


dh 

dy 


0 , 


f+C 

h 

dh 

^ t: (0-«/2) 


=const 


dh 


d \ x n^) 


( 12 ) 

(13) 


ZONAL SLOPE 



Figure 3. Diagram of Zonal Shear Scenarios 


4. Orientation Cases 

For the purposes of simulating initial flows for various orientations, we chose to 
retain the flow direction relative to the computational domain and re-orient the basin in 
order to simplify our analysis and allow us to induce periodic boundary conditions along 
the flow. To simplify our post-processing, we retain the flow direction and shape and re- 
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orient the basin. Future experiments may wish to try alternative methods of measuring 
and calculating orientation output across the flow for comparison. Rotating the basin 
involved converting from Cartesian grid, which we used on our other runs, to a spherical 
grid. The orientation was controlled by using Eulerian coordinates to turn the basin 
around the z-axis using the scale 1 degree latitude =111 km. This gave us a dx and dy of 
-1.8° vice 2 km. Figure 4 depicts the orientation scenarios and their respective angles. 



Figure 4. Diagram of Orientation Case Scenarios 

C. DIAGNOSTICS OF EDDY HEAT FLUX AND DIFFUSIVITY 

In order to evaluate the impact of topographic effects on mesoscale variability, we 
must anticipate what variables that our numeric model outputs that we would expect to 
see change. This section outlines some of the rationale behind our diagnostic tools based 
on the relationships to each other and to how we would expect them to be impacted under 
different conditions. Because of our model setup with a current that is strictly zonal 
(except for the “orientation” runs), meridional motion is mostly caused by the baroclinic 
instability. Therefore, our measurement of any transient displacement in y-direction will 
reflect the intensity of mesoscale variability. This is important because in order to 
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determine the effects of mesoscale variability, it should be clearly isolated from any 
direct effects of background flows. To accomplish this, we simplify the eddy equations of 
motion and make our initial background currents purely zonal flow. Inherent in this 
simplification is the lack of any rotational or non-meridional transport and so divergent 
and rotational components are not isolated (Radko et al. 2014). Note that for the 
orientation case, this simplification is only possible by applying values of x and y as 
“along-flow” and “across-flow,” respectively. Additionally, non-eddy effects are 
removed by removing mean flow from our v-velocity. 

As previously mentioned, mesoscale eddies are an important component in the 
transfer of heat from the equatorial region to the poles. This means eddies are also large 
contributors in maintaining the heat balance for our planet. Because the stratification in 
our idealized basin is maintained through temperature and not salinity, we are able to 
truncate our equation of state to essentially replace density gradient with thermal 
gradient. This allows us to approximate density with temperature, and therefore 
temperature becomes a key diagnostic variable. MITgcm allows you to output both 
temperature and heat flux in both x and y directions (Menemenlis and Fukumori, 2005). 
Because we are taking our values over time and we can average gridpoint values over the 
entire model area and over all the equilibrium time periods in order to get one value for 
the average temperature flux (Qj) over the period from when the model reaches baroclinic 
equilibrium to the last integration output. 

Q f = FT (14) 

If the heat flux is required, it can be trivially computed from the temperature flux 
as follows: 

Qhf = P c P Q f > (15) 

where C p «4000J kg~ l K~ l is the specific heat and p ~ 1 000kg m 3 is the density of 

seawater. Thermal diffusivity (K/j can be obtained by dividing the temperature flux by 
the thermal gradient, as indicated in Equation (16). 
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(16) 


dT 




dy 


K t 


T -T 

1 iH 1 Jj -1 


2 dy 

Qf_ 

T y 


Finding the overall thermal gradient in Equation (17) from the model is 
accomplished by averaging the temperature gradients for each gridpoint in space and 
time. The advantage behind diffusivity is that it can easily be compared with alternative 
estimates in literature and analytical models, like GRB. To maintain a statistically steady 
state, a relaxation condition is enforced on the surface and bottom in order to prevent 
unrealistic build-up of warm and cold temperatures. Additionally, all diagnostic variables 
are calculated only after removing the boundaries as far out as 60 km. This was necessary 
to be able to focus on what was occurring within the deep water basin and to avoid 
contamination from boundary layer dynamics. Think of this as removing the wax layer 
from around the cheese in order to prevent a nasty taste in your mouth. 
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III. RESULTS 


A. DIRECT VERSUS REMOTE SHEAR OVER MERIDIONAL SLOPE 

We compared the initial meridional slope runs conducted in the 1100-meter basin 
to determine any topographic effects on eddy heat flux. We also investigated the remote 
scenario, where the flow and shear were not in direct contact with the seafloor, for the 
effects of varying the slope. There is a remote effect on eddy heat flux even without 
direct interaction. However, the mean total heat flux is over three times greater in the 
direct case than in the remote case for an equivalent slope (s = 2.62xl0' 4 ) after the model 
reaches an equilibrium. In both cases, a south slope of the same absolute value as the 
north slope yields a greater meridional heat flux as seen in Figure 5. 



0.012 
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0.004 
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0.004 
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Heat Flux mean • Direct & Remote, North & Slope w/Cd 


500 1000 1500 

Days 


Figure 5. Direct Versus Remote Case Temperature Flux Comparison at 

Full North and Full South Gradients 


In our direct case, variations in the meridional bottom slope are expected to affect 

variation in total mean heat transport. Chen and Kamenkovich (2013) observed in a two- 

layer model that meridional slope changes the PV gradient causing a stabilizing effect. 

The results of the 1100-meter model indicate that in the direct shear case, the largest 

mean meridional heat flux and diffusivity occur for a flat seafloor (zero slope), which 
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would indicate that a bottom gradient of any kind reduces horizontal diffusivity (Figure 
6). However, calculating the meridional heat flux in the upper half of the basin, the 
remote case (Figure 7) yields the greatest flux for a slightly negative bottom slope. The 
same can be said for diffusivity (Figure 8) which shares the same curvature. 


Direct 1100m meridional slope scenario 



Figure 6. Direct Case Mean Meridional Temperature Flux and Diffusivity 

Versus Slope 


One of the explanations for the weaker total heat flux in the remote case could be 
that the remote case is more strongly affected by conservation of potential vorticity. 
Perhaps this is due to the effect of the bottom drag, which may have a stronger impact on 
direct flows because of the solid boundary directly below the shear region; the latter 
feature is absent in the remote case. In order to explain the pronounced asymmetry of 
temperature flux in the remote case (Figure 7), we now attempt to develop a simple 
analytical theory predicting the slope that results in maximal heat transport. This slope, 
which will be referred to as the optimal slope hereafter, is computed based on the analysis 
of the equation for conservation of potential vorticity (Q pv ), given by 


Q 


pv 


s + f 

h 


g + f 0 + Py 

h 0 -sv 


= const. 


(17) 
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Remote 1100m meridional slope scenario 


xIO' 3 Interpolated Temp Flux vs Slope - Upper half 



Slope xio' 4 


Figure 7. Eddy Thermal Flux in Upper Half of Remote 1100-meter Case 2 


Remote 1100m meridional slope scenario 
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Figure 8. Thermal Diffusivity in Upper Half of Remote 1100-meter Case 

Conservation of Q pv implies that the meridional displacement of water columns is 
adversely affected by both the |]-cffcct and by the slope individually. However, when 
both effects operate concurrently, they can counteract each other and thus have minimal 

2 MITgcm thermal flux output was plotted along with the calculated thermal flux for validation. 
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Thermal Diffusivity vs Slope - upper half 











impact on the meridional displacements of water columns. By assuming that this 
cancellation occurs when 


/ 0 + [iy oc \ - sy 

we can roughly identify an optimal slope ( S opt ) based on P, f, and the slope s, where 
/ = lxlO’ 11 m-s' 1 ,/) = lxlO" 4 s' 1 , and ho is the maximum depth of the water column. This 
assumption can be rationalized as follows. The conservation of Q pv states that if the 
height of the water column changes, rotational vorticity (Q resists change, and therefore, 
Q pv is conserved through a north or south shift, which compensates for changes in f. On 
the other hand, if/ changes, the water column will want to compensate by stretching or 
shrinking. Equation (18) then immediately implies 

PK 

(19) 


S op , = ■ 


fo 


This formula for optimal slope, Equation (19) shows a calculated S opt of 
approximately -l.lxlO' 4 for ho= 1100 m, which is comparable to our model results. Using 
a splines interpolation, the value of the maximum slope based on the compared model 
runs is -0.858xl0' 4 . 

Using this same equation, the 3300-meter remote case can also be analyzed. 
Because the theoretical S opt is a function of the base depth, it is easily estimated at three 
times the 1100-meter case, (S op ,(ha=33 00 m) = 3.3xl0' 4 ). However, our results are a bit 
more complicated than our simple theory. This meant that additional southern, or 
negative slope values had to be run for this case in order to find the slope great enough 
for our optimal case. However, even with these extreme slope simulations, the optimal 
slope for heat flux has not yet been attained. Figure 9 illustrates that the general shape of 
the slope to (9/ difference is maintained, the S 0/) , does not appear to be within the range of 
slopes modeled. In fact, it appears that the increased heat flux levels off below -4xl0‘ 4 
rather than tapering off sharply like they did in the 1100-meter case. 

It initially appears that the optimal slope theory tends to output too gentle a slope 
for the /?=3300 m case. However, taking a look at the interpolated diffusivity plot in 

Figure 10, the diffusivity does begin to drop off at a S int of -4.54xl0' 4 ’ which is close to 
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our calculated S op t■ A clear difference between the diffusivities and heat flux can be seen 
at larger slopes, whereas in the slopes calculated in the 1100-meter case, the Kt and Qf 
curves were nearly identical. 


Remote 3300m meridional slope scenario 



Figure 9. Eddy Thennal Flux in Upper Half of Remote 3300-meter Case 


Remote 3300m meridional slope scenario 
Interpolated Diffusivity vs Slope - Upper half 



Figure 10. Thermal Diffusivity with S opt and S int , Upper Half of Remote 

3300-meter Case 
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B. ZONAL SLOPE 


The zonal slope case was analyzed in the same way as the meridional slope case 
in order to compare the direct and remote scenarios using the two separate slopes. Results 
show that in both the direct and remote cases, and in both slopes, the net heat flux was 
equivalent to that of a flat bottom case. Singling out the £2/ observed on either the east or 
west side of the apex, we discover that they are noticeably large and erratic, which may 
be corrected if the turbulent fluxes were subtracted from the total flux. However, it is 
apparent to see that the southward mean flow of the western, or “uphill,” side of the basin 
is fully countered by the northward flow on the eastern “downhill” side. In Figure 11, the 
black and gold lines represent both the flat bottom case equivalent and the net Qj for the 
full water column for which the values are nearly identical. This was true both for the 
direct and remote cases at both slope angles where the calculated difference between the 
net zonal and flat bottom fluxes was negligible. 


Zonal Flux Comparison for Remote shear at h=0.5 



Figure 11. Zonal Slope Flux Comparison 


C. ORIENTATION 

As a reminder, diagnostics for the orientation runs were calculated in the frame of 
reference associated with the initially imposed basic flow. In particular, we look at the 
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cross flow heat flux, diffusivity, and the velocity and magnitude of the mean current. 
Figures 12 and 13 illustrate the different CFF and CF-DIFF patterns for each value of 0. 
As 0 increases, the CFF is dampened until (9=90° at which point it rises again. At 
(9=180° it reaches a maximum CFF, three times that of the zonal flow at 0=0° and in the 
polar opposite direction: southward flux instead of northward. Similar to larger slopes in 
the 3300-meter zonal flow model, the orientation models’ CF-DIFF and CFF plots do not 
have the same shape. CF-DIFF, for example, reaches its minimum at (9=0° and its 
maximum at (9=120°. This is 60° out of phase with the CFF curves of both values. The 
results suggest that diffusivity is offset from fluxes when oriented flow is involved. 
Additionally, the diffusivities are approximately twice the typical values observed in the 
ocean. Appendix A (Figure 19) depicts Figure 12 with orientation angles for better 
reference. 


Remote 3300m Orientation scenario 
xio" 3 Cross-Flow Temp Flux vs 0 - Upper half 
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0 degrees from X-axis 


Figure 12. Orientation Scenario CFF 
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Remote 3300m Orientation scenario 
Thermal Diffusivity vs 9 - Upper half 

12000 

* 

10000 
8000 
£ 6000 
4000 
2000 , 

0 

0 30 60 90 120 150 180 

9 degrees from X-axis 

Figure 13. Orientation Scenario CF-DIFF 

Mean velocities were also computed in order to identify potential changes in the 
properties of background flow. In Figure 14 (right), the angle of the mean flow is plotted 
in relative to the equator, or x-axis (magenta) and relative to the orientation of the initial 
flow (cyan). Not surprisingly, there is a clockwise veer for each of the non-zonal flows, 
which are all comparatively similar, roughly in the range of 4-12 degrees. Both of the 
zonal flows (0° and 180°) maintain their mean flow angle with the direction of 
background flow, ffowever, the magnitude differences between the mean flows at 
different values of 0 is noteworthy. Note in Figure 14 (left) that the mean flow velocity 
for 0=0° is the initial value of the mean current and is therefore approximated. Initially, 
the sharp decrease in mean flow velocity begins like the CFF, dropping significantly and 
rising again all the way up to 0=150°, it then suddenly drops again around 0=180°, 
unlike the CFF. 


Thermal Diffusivity vs 9 - Upper half 
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mean flow angle H 


Remote 3300m Orientation scenario 
Angle of mean flow - Upper half 
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Remote 3300m Orientation scenario 



6 degrees from X-axis 


Figure 14. Angle of Mean Flow (right) and Magnitude of Mean Flow (left) 
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IV. COMPARISON WITH GRB 


This chapter focuses on comparing our 3-D model data to GRB. First by briefly 
explaining the theory of GRB, then by comparing our 3-D MITgcm results to the GRB 
results for our 3300-meter remote case. 

A. THEORY 

Oceanographers have often attempted to explain the non-linear effects of BI using 
analytical models. GRB attempts to explain BI by linking the growth rates of the primary (Xi) 
and of the secondary (a 2 ) modes of instability. Such analysis can be insightful because, as the 
primary mode increases in intensity, the secondary mode, which is a function of the primary, 
intensifies and eventually dampens the primary mode. In this theory, the two modes 
“compete” and, in doing so, an equilibrium state is created where the instabilities finally 
match each other. This “balance” can be calibrated using an empirical constant, (C). 


(a) 

North 


(b) 

North 



Figure 15. Schematic of GRB Model. Source: Radko et al. 2014. 
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One of the features of this theoretical model is that the primary modes are a 
function of vertical shear, bottom drag, stratification, and [3-effect, whereas the secondary 
mode is a function of all of those plus the primary mode amplitude. The GRB equation 
can be summed up by the following: 

4 = C\ (20) 

In theory, the primary mode would grow continuously if not for the secondary 
mode hindering its development. Likewise, if the secondary instabilities (k 2 ) were to 
dominate, the primary instabilities would be overly-dampened which would ultimately 
result in no eddy formation. This means the two modes must reach a balance at some 
point, at which the output diagnostics can be time averaged to see the mean results. 
Examples of this instability and its progression can be seen in a horizontal section of the 
numeric models as shown in Figure 16. The Growth Rate Balance theory is discussed in 
greater detail by Radko et al. (2014). For this paper, numerical results from 3-D MITgcm 
model runs will be used to compare to a 2-D Growth Rate Balance analysis. Figure 16 
includes three time steps of potential vorticity from a numeric model that illustrates the 
progression of initial growth rate (left), the development of secondary instabilities 
generating wave-like patterns (center), and the resulting transient irregularities (right). 




Figure 16. Three Timesteps of Model Potential Vorticity. Adapted from 

Radko et al. 2014. 
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B. RESULTS 


In determining the best fit for the GRB model, Radko et al. (2014) investigated 
several values of the non-dimensional empirical constant, C, that would best fit the 
Phillips (1951) and Eady (1949) models of linear instability. A range of 3.5-4 was 
determined to be the best fit, which was calibrated using PV flux. The 3-D models to 
which they compared to had spatial scales similar to our model, lacked sloping 
topography and the shear was in the direct contact with the seafloor. For this thesis, the 
3300-meter remote case was plotted against calibrated values of PV diffusivity 
determined from GRB models for each slope at values of C from 2 to 5. Interestingly, the 
GRB models predict a slightly northern optimal slope. Only the diffusivities calculated 
from our flat bottom (zero slope) 3-D MITgcm model compared well to the GRB model, 
which was just under the C=3.5 plot. As seen by the 3D versus C values in Figure 17, 
within a range of slopes from +2.6xl0' 4 to -1.4xl0" 4 the various diffusivities did stay 
within the boundaries of 2 <C <_5, but for any slope outside of that range, the analytic 
model was not ideal. 



Figure 17. 3-D MITgcm 3300 m Diffusivities compared to GRB 
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Ultimately, the gradual bell curve of the GRB runs does not match the curve of 
the 3-D MITgcm plots (Figure 17). If we look at a value of C = 3.5, which is consistent 
for our flat bottom case, the predicted diffusivities from the GRB theory are too low for 
negative slopes, and it appears to be too symmetric along the y-axis. It appears that, 
although the GRB model compared well to the Eady model in Radko et al. (2014), where 
C = 3 fit for each run, for our model, changing the slope of the sea floor creates 
variability that the GRB model does not predict. Comparisons with the quasi-geostrophic 
model can be seen in Appendix A (Figure 20). 

Several factors may be able to explain why GRB does not do well with 
topographic slope. One may be the simulation depth. The 3-D model in Radko et al. 
(2014) was 1100-meters deep. Here, we are comparing the 3300-meter model. Second, 
the assumptions in the GRB model may not accurately account for topography in the 
basic equations. Although the GRB model accounts for topographic variation, it may not 
be able to represent more complex geophysical domains. Third, there may be a factor in 
developing BI that GRB is not able to account for. Ultimately, however, this is a good 
indication of why we need a better understanding of the effect bottom topography has on 
dynamics. 
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V. DISCUSSION AND CONCLUSIONS 


The first significant conclusion or this research is that the remote case shear 
scenario acts differently than the direct case with regard to the effect on BI. First, the 
intensity of heat transport in the direct case is greater given the same shear and slope. 
This is to be expected since the shear force covers twice the amount of volume as in the 
remote case. What is unexpected is that the meridional slope acts on the direct shear case 
differently than on the remote case. Any meridional slope dampens the heat flux in the 
direct case, whereas in the remote case there is a non-zero optimal slope for maximum Qf 
It should be noted that this is a large-scale experiment, but similar effects may also occur 
in smaller scale scenarios. This could be important in straits where there is a meridional 
slope with seasonally changing shear depths: a scenario where the difference in heat 
transfer may prove more complicated than typical parameterizations account for (Spall 
and Chapman 1998). 

The distinct optimal slope of the remote case is perplexing since, based on P- 
effect alone, the direct case should also have a non-zero optimal slope. One possibility is 
that in the direct case, shear acting directly on the motionless seafloor is more influenced 
by bottom drag. Since the shear is directly impacted by the bottom drag, the frictional 
forces due to any increased slope overcome the variation caused by the P-effect, therefore 
making the flat bottom the most efficient for Qf On the other hand, the shear in the 
remote case is shielded from the bottom by an initially quiescent layer, which means the 
effect of bottom friction is weaker. This could prove useful in determining 
parameterizations of flux and diffusivity based on the depth of the current 

With increased slope, the diffusivity curve and heat flux curves no longer match. 
This is illustrated quite clearly in the 3300-meter Kj and Qf where the interpolated 
optimal slope for K T is of smaller value than that for Qf. Similarly, the orientation cases 
show a maximum angle for Qf that is not the same orientation angle as those with 
maximum Kt. In cases of extreme slope or orientation, heat flux and diffusivities are no 
longer proportional. This result is significant in that it justifies making model 
comparisons using multiple diagnostics. 
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In the orientation cases, the thermal gradient is shifted along with the flow to 
where it runs east—west and south—north with the same applied Coriolis force, 
effectively changing the realm of what we see on earth. With the Coriolis force being 
weaker at the equator, one might expect that the reason for this intense (Qf in our (9=180° 
model is due to our horizontal thermal gradient. If we look at this in terms of what would 
happen if warmer water were to the north vice the equator, and neglected any changes in 
atmospheric dynamics that would occur from this, we would have a thermal/density 
gradient and Coriolis gradient running parallel. Unlike in nature where the wannest water 
temperatures occur where there is the least Coriolis force, orientations above 120° have a 
stronger Coriolis is in the same region as the wanner water. Perhaps this scenario 
generates more eddy motion / mixing in the north that what we observe near the equator. 
This may act as a “heat engine” that, along with thennal wind, promotes larger amounts 
of CFF. Perhaps in the future research is needed to validate this concept. 

For the GRB model, the results illustrate that the assumptions built into this quasi- 
geostrophic formulation of the GRB model only capture very approximate qualitative 
behavior. The 3300 m 3-D model optimal slope falls in line with comparisons done 
between the basic two-layer QG model, as seen in Appendix A (Figure 20). However, the 
QG slope and magnitude are still half that of the 3-D scenario. This reaffirms that bottom 
slope does affect BI and that two-layer and analytic models may lack the ability to 
properly account for meridional-sloped topography. 
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VI. RECOMMENDATIONS 


There is a broad range of possible studies that can be done to continue this 
research as many other practical comparisons could be made. For example, changing the 
strength of shear, or comparing heat flux and diffusivities at different slopes to see when 
they match and when they begin to diverge. Below is a short summary of topics that one 
may wish to explore, but there are many other possibilities that would undoubtedly be 
directly relevant to oceanography and geophysical fluid dynamics. 

1. The effects of mesoscale and sub-mesoscale bottom features on meridional 
transport: This study could use the same numeric model setup but create one or many 
bottom features replicating underwater mountains or atolls with scales from 10 to 100 km 
and heights extending through a significant fraction of water depth. A comparison could 
then be made between the total meridional heat and PV flux with and without these 
features and according to their size and comparing them to the uniform slope data in this 
paper. 

2. Examining the effects of different shear patterns to determine what shear 
characteristic plays the dominant role in the dynamics of Bl. This can be accomplished 
by varying the speed of the surface current and thus changing the gradient of the flow 
over depth, and by introducing various shears in different water depths. 

3. Continued work on the orientation cases to improve our understanding of 
physical processes controlling the variability of heat transport with the flow direction. 

4. Detennining the effect of eddies in remote and direct case on sea surface height 
and temperature to determine and correct systematic errors in satellite bathymetry 
calculations. This experiment can be initiated by attempting to recreate a height-varying 
topography on the model seafloor from the sea surface data of the ocean using the model 
output data post-equilibrium. 
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APPENDIX A. COMPLEMENTARY FIGURES 
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Figure 18. 3-D Captures of 1100-meter Direct and Remote Model Runs 
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Meridional Diffusivity (m 2 /s) 


Remote 3300m Orientation scenario 
xio 3 Cross-Flow Temp Flux vs 0 - Upper half 
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9 degrees from X-axis 


Figure 19. Orientation CFF with Orientation Diagram for Visual Reference 
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Figure 20. GRB with 3-D 3300-meter and Quasi-Geostrophic Comparison 
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APPENDIX B. DATA TABLE 


Below is the data collected from various model runs. Orient_ represents the 
orientation runs, 3D_ represents the 1100 m zonal flow remote case, and 8D_ represents 
the 3300 m zonal flow remote case. 


JOB 

Depth 

slope 

theta 

Qf 

Qf MITgcm 

KT 

orient30N 

3300 

0 

30 

6.164E-03 

6.155E-03 

1.740E+03 

orient60N 

3300 

0 

60 

3.847E-03 

1.197E-02 

2.074E+03 

orient90N 

3300 

0 

90 

2.794E-03 

1.897E-02 

3.549E+03 

orientl20N 

3300 

0 

120 

3.424E-03 

1.960E-02 

9.046E+03 

orientl50N 

3300 

0 

150 

5.827E-03 

1.897E-02 

1.155E+04 

orientl80 

3300 

0 

180 

9.665E-03 

1.611E-02 

1.140E+04 

3D SN RL 5 

1100 

-2.600E-04 

0 

6.637E-03 

6.715E-03 

1.762E+03 

3D SN RL 25 

1100 

-1.300E-04 

0 

7.486E-03 

7.440E-03 

2.135E+03 

3D SN RL 125 

1100 

-6.550E-05 

0 

8.104E-03 

8.030E-03 

2.282E+03 

3D RL FB 

1100 

0.000E+00 

0 

5.422E-03 

5.516E-03 

1.468E+03 

3D NS RL 125 

1100 

6.550E-05 

0 

4.231E-03 

4.185E-03 

1.065E+03 

3D NS RL 25 

1100 

1.300E-04 

0 

3.148E-03 

3.110E-03 

7.610E+02 

3D NS RL 5 

1100 

2.600E-04 

0 

2.059E-03 

1.998E-03 

4.766E+02 

8D SN RL100 

3300 

-5.244E-04 

0 

1.075E-02 

1.078E-02 

3.507E+03 

8D SN RL 75 

3300 

-3.933E-04 

0 

1.061E-02 

1.056E-02 

3.521E+03 

8D SN RL 5 

3300 

-2.600E-04 

0 

1.023E-02 

1.013E-02 

3.328E+03 

8D SN RL 25 

3300 

-1.300E-04 

0 

9.067E-03 

9.148E-03 

2.791E+03 

8 SN RL 125 

3300 

-6.550E-05 

0 

7.438E-03 

7.340E-03 

2.194E+03 

8D RL FB 

3300 

0.000E+00 

0 

6.164E-03 

6.155E-03 

1.740E+03 

8D NS RL 125 

3300 

6.550E-05 

0 

5.483E-03 

5.462E-03 

1.491E+03 

8D NS RL 25 

3300 

1.300E-04 

0 

5.347E-03 

5.316E-03 

1.435E+03 

8D_NS_RL_5 

3300 

2.600E-04 

0 

4.091E-03 

4.059E-03 

1.050E+03 
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